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Abstract 

On the basis of physical considerations we propose a one-dimensional 
discrete lattice model for the density relaxation of granular materials under 
tapping. Solving the difference equation numerically, we find a logarithmic 
time-dependence of the density relaxation. This is in agreement with exper- 
imental results of Knight et al. [Phys. Rev. E 51, 3957(1995)]. The origin of 
this anomalous relaxation is elucidated analytically by solving the equation 
of its continuum version asymptotically in time. 
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I. INTRODUCTION 



The slow relaxation of density in granular materials under tapping has attracted much 
theoretical interest recently since Knight et al. carried out systematical experimen- 
tal measurements of density as a function of time for a vibrated granular material. The 
experimental data of Knight et al. can be most satisfactorily fitted using a functional 
form: 

where p(t) is the average density in the system at time instant t. Here poo, ^Poo, G and r are 
constant parameters. This four-parameter fit was obtained from experimental data without 
theoretical motivation. Previous theoretical models predicted different functional forms. In 
the model of Barker and Mehta [0], particles can relax both independently, as individual 
particles, and collectively, as clusters. Their model leads to a sum of two exponential terms, 
which is not in accord with the logarithmic form of Eq.(l). The model proposed by Hong 
et al. p is based on a diffusing void picture. It predicts a power-law dependence of height 
reduction as a function of time, 6h ~ t^ with z = 1, which means that the density increases 
linearly until saturation. 

Using frustrated Ising models, Herrmann and his collaborators gave a numerical 
confirmation of the logarithmic density relaxation as Eq. ([l|) in their Monte-Carlo sim- 
ulations. In their models, frustration plays a crucial role. They claim that frustration is 
generated in granular materials by the steric constraints imposed by the hard-core repulsion 
of neighboring grains and the subsequent interlocking. The link of the frustrated lattice gas 
models to real granular packing is explained in Refs. ffj-Q. 

In this paper, we propose a one-dimensional discrete lattice model for granular com- 
paction under tapping based on a simple physical picture. Numerical solutions to the model 
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equation indicate that the density relaxes in a logarithmic manner, precisely as Eq. (1). We 
also solve the continuum version of the model analytically in the asympototic limit {t — > oo) 
to confirm Eq. (1). 

This paper is organized as follows. Section II presents the model equation. Numerical 
and analytical solutions to the equation are given in Sections III and IV respectively. Section 
V is devoted to discussion. 



II. MODEL EQUATION FOR DENSITY RELAXATION 

The density of particles p(r, t) should satisfy the continuity equation: 

^p(r,i)+V-J = (2) 

where J is the current density. We may rule out any simple density-gradient term in the 
current density which leads to isotropic diffusion, since the thermal energy is too small to 
triger the motion of grains. The current density may be written as a product of density and 
velocity: 

J = p(r,i)v (3) 

To specify the the functional form of velocity, let us consider a one-dimensional lattice model 
along z-direction. Gravity is along the negative z-axis. In the absence of tapping, motion of 
grains under the action of gravity is possible only when the geometric constraint is satisfied, 
i.e., enough void space below the grain. If p{z,t) is the particle density, the void density at 
the site just below is 1 — p{z — 1, t). Therefore, the velocity is nonzero only when the ratio 
a — i_p^^^'% less than 1. This means a step function for the velocity v: 



V ^ < 



—vq , a < 1 

(4) 

, a > 1 



in the absence of tapping. The effect of tapping in granular compaction is to overcome the 
bottleneck effect, i.e., to make the geometric constraints not as strict as the above. Motion 



of grains is also possible when a > 1 with the help of tapping. We may therefore replace 
the above step function with a peak function of an exponential form: 

V = -De-"/"' (5) 

where 7 and D are constants dependent of the vibrational intensity F = Au'^/g with the 
amplitude A and the frequency u of oscillation and g the gravitational constant. The precise 
F-dependences of D and 7 are beyond the scope of the present theory. Substituting Eq. (H) 
and Eq. (H) into Eq. (||), we obtain our model equation for density relaxation: 

|p(^,^) = ^^^(^,^) (6) 

with 

where ^ is infinitesimal in the continuum limit. This model is defined more definitely on a 
lattice: 

^P^ = D{Fiz + l,t)-Fiz,t)) (8) 
where F{z,t) is given by Eq. (0) putting ^ = 1. 



III. NUMERICAL SOLUTIONS 

In order to solve numerically Eq. we supply two fixed boundary conditions, p(0, t) = 
1 at the bottom and p{L + l,t) = at the top of the system. The following discrete version 
of Eq. (ID ensures no flux at the boundaries, i.e., conservation of particle density as a whole 
is ensured at each time step 

p{z, t + l)= p{z, t) + D{F{z + 1, t) - F{z, t)) (9) 

with F{z,t) given by Eq. (|^) with ^ = 1. 
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We update {p{z, t), z = 1,2, ■ ■ ■ , L} according to Eq. @ from random initial conditions 
with {p{z,0),z = 1,2, ■ ■ ■ , L} equal to random numbers distributed uniformly from pi to 
Pu- Fig. 1 shows snapshots of density profiles at three different time steps. We see that 
as particles move downward voids move upward and pile at the top. As soon as we start 
updating from a random initial configuration, a sharp interface between a particle phase and 
a void phase begins to emerge. Its position can be easily identified as the location where 
density changes sharply from non-zero value to zero. Since the position of the interface takes 
only integer values on the lattice and information extracted from it is therefore limited, we 
resort to calculate the position of center of mass 

H{t) = j\p{z,t)dz 

= jZzp{z,t). (10) 

z=l 

We take it for granted that the average density measured in experiments is proportional 
to the inverse of H{t). We determine the proportional prefactor as follows. The minimum 
value of H can be obtained from a density distribution of a step function, i.e.. 



p{z,t) = l ' - (11) 




which is actually a static solution to the model equation Eq. (^. Here B is determined by 
the density conservation and is a time-independent quantity: B = p{z,t)dz. The value 
of H corresponding to Eq. (|ll|) is \B^. We therefore use the following expression for the 
average density 

Pit) = (12) 

Starting from random initial configurations, we discard the data for initial transient time 
period and do not make computation of p{t) until at the interface between the particle phase 
and the void phase the density jumps from zero to some reasonably finite value. Our density 
of unity corresponds to the densest close packing and it sets the scale for density. For the 



case of sphearical particles in reality, the densest close packing has a density of about 0.64 
while the mechanically least stable configurations have the most loose value of about 0.55 



model. The time step when we start calculating p{t) in Fig. 2 is set to be the time origin. 
It is the instant when the density gap at the interface is about 0.8 in our density scale. The 
numerical data in Fig. 2 can be very well fitted by the logarithmic form of Eq. (|l]), which 
is also shown in the figure. 

Fig. 3 displays two density profiles at different time steps for the particle phase while 
zero density in the void phase is not plotted. We find that the density profile can be fitted 
by the following expression 



where pi{t), s{t) and f3{t) are time-dependent parameters. The fit to the numerical data is 
also plotted in Fig. 3. In the whole range the density is step-function-like 



where b{t) is the location of the interface. 

Eq. (p^ ) is different from the generalized Fermi-Dirac distribution proposed by Hayakawa 
and Hong [|1^] and found by Herrmann and his collaborators in simulations of frustrated Ising 
models ^j^. The generalized Fermi-Dirac distribution does not have a singularity in the 
density profile but changes continuously from larger values to zero rapidly within a boundary 
layer whose width is determined by one of the parameters in the distribution. On the other 
hand, our model gives a sharp interface between the particle phase and the void phase and 
therefore the model does not permit any density value smaller than the most loose density 
for granular compaction. In this sense, we may say that in addition to the granular solid 
phase and the void phase the generalized Fermi-Dirac distribution allows a granular gas 
phase where density changes drastically |]10|,^Q. Our model describes only the granular 



(in our scale 0.859). Fig. 2 shows one typical plot of p(t) calculated from the lattice 



(13) 




(14) 
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compaction under tapping where the grains do not fly as in a gas (and therefore no inertia 
terms appear in the velocity expression Eq. (^). 



IV. ANALYTICAL SOLUTION IN THE ASYMPTOTIC LIMIT 

As we observed in the simulations, there is always a sharp interface between the particle 
phase and the void phase. We may therefore generally postulate the solution to Eq. (||) as 

p{z,t) = a{z,t)e{b{t) - z) (15) 

where 0(x) is the step function such that 0(x) = for x < and 0(x) = 1 for x > 0. 
We now determine the time dependence of the location of the interface b(t) by solving the 
equation in its continuum version, Eq. (||). Substituting Eq. (p!5| ) into the left hand side of 
Eq. (P), we obtain 

^ = d{z, t)Q{b{t) -z) + a{z, t)h{t)5{h{t) - z) (16) 
ot 



where dot stands for time deriative. The G function in Eq. (|T^) gives us the following 
F{z,t) inEq. (0) 

F{z, t) = a{z, t)e{b{t) - ^)e-'^(-'*)/(^(i-'^- (^'*))) (17) 

= f{a{z,t),a_{z,t))e{b{t)-z) (18) 

where a-{z,t) stands for the value of a just below z, i.e. a{z — ^,t) with infinitesimal ^. 
Here 

f{a{z,t),a4z,t)) = a(2^t)e-"(^'*)/(^(i-"-(^'*») (19) 

Thus, the right hand side of Eq. (^) reads 

D^^^ = D^^m -z)- DfSm - z) (20) 

Integrating the right hand sides of Eq. (jl^) and of Eq. ( PO] ) from z = b{t)—A to z = b{t)+A 
and taking the limit of A —> 0, we obtain 



a{z,t)b{t)U=m = -I^a(z,t)e-'^(^'*)/(^(i-'^-(^'*«) (21) 

Now we make an approximation. From the simulations we know that a{z, t) depends on 
z very weakly. The parameter (3{t) in Eq. (|13D decreases as time increases. When t oo, 
(3^0. We therefore make the following aussmption in the asymptotic limit of t — » oo 

a{z,t) = a{t). (22) 

The conservation of total density is then expressed as 

a{t)b{t) = B (23) 

with the same B as in Eq. (|ll|)- 

Using Eq. (|l|), Eq. (|2|) and Eq. (|2|) we have 

Bd = Da^e'^^) (24) 

Letting a = 1 — e leads in the asymptotoc limit (e 0) 

e = — ce~^ (25) 



where c = -^^^^ is constant. Putting x = 1/e, we have 



X = cx'e"^'/^. (26) 

Now we integrate Eq. (|26D for t from to t and for x from to where po and p 
are the average density (i.e., a(t)) of the particle phase at time and t respectively 



x-^e^'/^dx = cj dt. (27) 

The asymptotic solution of Eq. (pT]) is readily obtained by using the formular 

POO p—y 

Eii-x) = - / dy (28) 

Jx V 



and the expansion for x >> 1 
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E^{-x)=e-[-l + -^+0{l^]. (29) 



The final form for Eq. is given by 

7(1 - p)2ei/(7(i-p)) =ct + d (30) 
with d = 7(1 — p(jYe^^^'^^^~''°^\ This can be rewritten as 

In 7 + 2 ln(l - p) + — ^ — - = ln(ct + d) (31) 
7(1 - p) 

where, however, the first two terms are neghgible compared with the third one for the left 
side, so that we obtain the final form of the average density 



This is precisely the same form as Eq. (||). The form of Eq. (|32D perfectly agrees with the 
simulations as shown in Fig. 2. 



V. DISCUSSION 

We have shown that the discrete lattice model Eq. (H) exhibits an anomalously slow 
density relaxation. Although our model is quite simple, it contains essential features of 
granular materials that granular particle can move only when there is enough void below 
the particle. Tapping the system weakens this constraint and we have taken into account this 
property by assuming the velocity as Eq. (^. In this way, the logarithmic time-dependence 
has been obtained, which is consistent with experiments. To our knowledge, no model 
equation in terms of space-time coordinates, which produces a logarithmic relaxation, has 
been available so far in granular systems. 

The analytical study indicates that the logarithmic relaxation originates from the factor 
exp[— p(2;)/(l — p{z — 1)] in the velocity Eq. (^. When p{z) ^ p{z — 1) ^ 1, the velocity 
becomes extremely small whereas it is finite when p{z) ^ 0. Thus we expect that the 
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logarithmic relaxation occurs as far as the system processes this property independent of the 
detailed form of the velocity. Actually, we have confirmed by simulations that the velocity 
such as V = D/{1 + e("-^)/^) leads to the essentially same relaxation and the existence of a 
sharp interface in the density profile. This fact implies a universal feature of the logarithmic 
relaxation not only restricted to granular materials but also for other systems where velocity 
vanishes abruptly like an essential singularity when the density exceeds some threshold value. 
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FIGURES 

FIG. 1. Snapshots of density distribution for three different time steps for system with 
L = 1024, 7 = 1 and D = 1. (a) At t = 0, the density is randomly distributed witli a mean 
value of 0.50. (b) Att = 1136 when the density jump is about 0.8 at the interface between the par- 
ticle phase and the void phase, (c) At t = 1024, 000 density distribution is very weakly dependent 
of position in the particle phase. 

FIG. 2. Time dependence of the average density. Here t = is shifted to the time instant when 
the calculation of p{t) starts. Data points are numerical results from the 1-D lattice model with 
L = 1024, 7 = 1 and D = 1. Best fit using Eq. (1) is also plotted (curve) with fitting parameters 
Poo = 0.989749, (5/9oo = 0.150325, G = 0.230741, r = 806.443. The agreement is remarkable. 

FIG. 3. Density profiles (data points) at two different time steps t = 10240 (lower) (start- 
ing from a random initial configuration) and t = 102400 (upper) obtained from the 1-D lat- 
tice model with L = 1024, 7 = 1.8, D = 1. Best fits using Eq. (13) are also shown (curves) 
with fitting parameters pi = 0.960768, s = 5.35083 x 10"^, /3 = 2.15145 for t = 10240 and 
pi = 0.967761,5 = 5.99141 x lO'^,/? = 1.8166 for t = 102400. 
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Fig. la -Peng and Ohta 




Fig. lb -Peng and Ohta 
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Fig. Ic -Peng and Ohta 
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Fig. 2 -Peng and Ohta 



